Creating maps with R is usually an easy task. However, there are cases where it can be hard to get a good representation, specially when involving maps where the internantional date line is centered on the map or projections involving the poles.

The main problem lies in the fact that different mapping file providers (specifically vector files aka “shapefiles”) uses different approaches, such as “breaking” the shapes on that coordinates (e.g cutting the Chukchi Peninsula) or just not providing data.

In this post I am going to fix the geometry of the shapefile provided by GISCO for the Antarctica and make some cool maps with the final geometry on orthographic projection. Let’s start!

# Libraries
library(tidyverse)
library(sf)
library(giscoR)
library(ggrepel)
library(rmapshaper)

Fixing the geometry

First, we get the map from GISCO:

antarct <- gisco_get_countries(year = 2024, resolution = 1, country = "ATA") %>%
  select(NAME = NAME_ENGL) |>
  # Ortho proj centered in the South Pole
  st_transform(crs = "+proj=ortho +lat_0=-90 +lon_0=0")

ggplot(antarct) +
  geom_sf()

And here we have it. This shapefile is provided with some lollypop-like cut that looks really ugly on a map. We are going to correct it manually with a few data wrangling. Basically we are going to:

  1. Identify and extract the polygon that represent “mainland” Antarctica.
  2. Convert it to coordinates (points).
  3. Identify those points that we want to remove and,
  4. Re-convert the coordinates to the polygon and replace the broken geometry with this fixed version.

We can convert from polygon to points as;

# Identify the max
ant_explode <- antarct |>
  st_cast("POLYGON")

nrow(ant_explode)
#> [1] 778

# Max polygon

ant_max <- ant_explode[which.max(st_area(ant_explode)), ]

coords <- st_coordinates(ant_max) |>
  as_tibble() |>
  # Add id for points
  mutate(np = row_number())


ggplot(coords, aes(X, Y)) +
  geom_point(size = 0.05) +
  geom_text(aes(label = np), check_overlap = TRUE) +
  coord_equal()

Here we can see that the offending points lies in the range of (8200, 9200). Let’s take a closer look:

test <- coords |>
  filter(np %in% seq(8200, 9200))

test |>
  ggplot(aes(X, Y)) +
  geom_point(size = 0.05) +
  geom_text(aes(label = np), check_overlap = TRUE)

Now we only have to play with the id (np variable) to capture the points we want to remove

Note that this cleansing has to be done individually for each shapefile. It is fairly simple but the code used here is specific for this shapefile.

# Final solution after some iterations...

test |>
  filter(np %in% seq(8289, 9130)) |>
  ggplot(aes(X, Y)) +
  geom_point() +
  labs(title = "To remove")


test |>
  filter(!np %in% seq(8289, 9130)) |>
  ggplot(aes(X, Y)) +
  geom_point() +
  labs(title = "To keep")

And now we can regenerate our shapefile with the right points:

# From coordinates to polygon
newpol <- coords |>
  as.data.frame() |>
  filter(!np %in% seq(8289, 9130)) |> # Removing offending points
  select(X, Y) |>
  as.matrix() |>
  list() |>
  st_polygon() |>
  st_sfc() |>
  st_set_crs(st_crs(ant_max))

ggplot(newpol) +
  geom_sf()


ant_max_fixed <- st_sf(st_drop_geometry(ant_max), geometry = newpol)

# Regenerate initial shape
antarctica_fixed <- bind_rows(
  ant_max_fixed,
  ant_explode[-which.max(st_area(ant_explode)), ]
) |>
  group_by(NAME) |>
  summarise(m = 1) |>
  select(-m) |>
  st_make_valid()

antarctica_fixed
#> Simple feature collection with 1 feature and 1 field
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -2576235 ymin: -2451273 xmax: 2683442 ymax: 2226821
#> Projected CRS: +proj=ortho +lat_0=-90 +lon_0=0
#> # A tibble: 1 × 2
#>   NAME                                                                  geometry
#> * <chr>                                                       <MULTIPOLYGON [m]>
#> 1 Antarctica (((-3390.457 -2450853, -3415.175 -2451201, -3594.403 -2451273, -36…

ggplot(antarctica_fixed) +
  geom_sf()

Plotting examples

With the right shape we can start mapping as we like. In this section I would try to reproduce some proposed designs of the Antarctic flag.

Graham Bartram’s proposal (1996)

Faily simple:

bbox <- st_bbox(antarctica_fixed) # For limits on the panel

antarctica_fixed |>
  ggplot() +
  geom_sf(fill = "white", color = NA) +
  theme(
    panel.background = element_rect(fill = "#009fdc"),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  ) +
  labs(title = "Graham Bartram's proposal") +
  coord_sf(
    xlim = c(bbox[c(1, 3)]) * 1.8,
    ylim = c(bbox[c(2, 4)]) * 1.4
  )

Emblem of the Antarctic Treaty

This one is a little more complicated, as we need the graticules forming a bullseye on the Antarctica. There are some additional challenges:

# Need graticules
grats <- giscoR::gisco_get_countries() |>
  st_transform(st_crs(antarctica_fixed)) |>
  # Specify the cuts of the graticules
  st_graticule(
    lat = c(-80, -70, -60),
    lon = seq(-180, 180, 30),
    ndiscr = 10000,
    margin = 0.000001
  )


ggplot(grats) +
  geom_sf()

We are going to merge the meridians so the South Pole is filled, as st_graticules() left a hole on it:

# Merge meridians
merid <- lapply(seq(-180, 0, 30), function(x) {
  df <- grats |>
    filter(type == "E") |>
    filter(degree %in% c(x, x + 180))

  df2 <- df |>
    st_geometry() |>
    st_cast("MULTIPOINT") |>
    st_union() |>
    st_cast("LINESTRING")

  sf_x <- st_sf(
    degree = x,
    type = "E",
    geometry = df2
  )
}) |> bind_rows()


grats_end <- merid |>
  bind_rows(grats |>
    filter(type != "E"))

And finally:

# Cut since some grats should be colored differently

antarctica_simp <- rmapshaper::ms_simplify(antarctica_fixed, keep = 0.005)
plot(antarctica_simp)

grats_yes <- st_intersection(grats_end, antarctica_simp)
grats_no <- st_difference(grats_end, antarctica_simp)

antarctica_simp |>
  ggplot() +
  geom_sf(fill = "white", color = NA) +
  theme(
    panel.background = element_rect(fill = "#072b5f"),
    panel.grid = element_blank(),
    axis.text = element_blank(),
    axis.ticks = element_blank()
  ) +
  geom_sf(data = grats_yes, color = "#072b5f", linewidth = 1) +
  geom_sf(data = grats_no, color = "white", linewidth = 1) +
  coord_sf(
    xlim = c(bbox[c(1, 3)]) * 1.8,
    ylim = c(bbox[c(2, 4)]) * 1.4
  ) +
  labs(title = "Emblem of the Antarctic Treaty")